Local structure of liquid carbon controls diamond nucleation 
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Diamonds melt at temperatures above 4000 K. There are no measurements of the steady-state 
rate of the reverse process: diamond nucleation from the melt, because experiments are difficult at 
these extreme temperatures and pressures. Using numerical simulations, we estimate the diamond 
nucleation rate and find that it increases by many orders of magnitude when the pressure is increased 
at constant supersaturation. The reason is that an increase in pressure changes the local coordination 
of carbon atoms from three-fold to four-fold. It turns out to be much easier to nucleate diamond in a 
four-fold coordinated liquid than in a liquid with three-fold coordination, because in the latter case 
the free-energy cost to create a diamond-liquid interface is higher. We speculate that this mechanism 
for nucleation control is relevant for crystallization in many network-forming liquids. On the basis 
of our calculations, we conclude that homogeneous diamond nucleation is likely in carbon-rich stars 
and unlikely in gaseous planets. 

PACS numbers: 



Most liquids can be cooled considerably below their 
equilibrium freezing point before crystals start to form 
spontaneously in the bulk. This is caused by the fact 
that microscopic crystallites are thermodynamically less 
stable than the bulk solid. Spontaneous crystal growth 
can only proceed when, due to some rare fluctuation, one 
or more micro-crystallites exceed a critical size (the "crit- 
ical nucleus"). An estimate of the rate at which critical 
nuclei form in a bulk liquid can be obtained from Classi- 
cal Nucleation Theory (CNT) 0. This theory relates R, 
the number of crystal nuclei that form per second per cu- 
bic meter, to AG cr j t , the height of the free-energy barrier 
that has to be crossed to nucleate a crystal: 
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Here k is a kinetic prefactor, T is the absolute tempera- 
ture and ks is Boltzmann's constant. The nucleation rate 
depends strongly on the height of the nucleation barrier. 
CNT predicts the following expression for the height of 
the nucleation barrier: 
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where 7^,5 is the liquid-solid surface free energy per unit 
area, Afj, is the difference in chemical potential between 
the solid and the supercooled liquid, and ps is the number 
density of the crystalline phase. The factor c depends on 
the shape of the nucleus, e.g. c = 167r/3 for a spherical 
nucleus. As the nucleation rate depends exponentially 
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on AG cr it, a doubling of jls may change the nucleation 
rate by many orders of magnitude. In general, the kinetic 
prefactor k in Eq. [1] can be estimated quite well Q • 

Because of the extreme conditions under which ho- 
mogeneous diamond nucleation takes place, there have 
been no quantitative experimental studies to determine 
its rate. Moreover, there exist no numerical estimates of 
A/u and 7^3 for diamond in supercooled liquid carbon. 
Hence, it was thus far impossible to make even an order- 
of-magnitudc estimate of the rate of diamond nucleation. 

In this Letter, we calculate the diamond nucleation 
rate R in liquid carbon at two state points {P= 85 GPa, 
T= 5000 K} and {P= 30 GPa, T= 3750 K} (points A 
and B in the carbon phase diagram shown in Fig. 1). At 
both state points, the liquid is supercooled by (T m — 
T)/T m w 25 % below the melting curve of diamond, 
with T m the melting temperature and = 6600 K and 
T r f t = 5000 K, respectively. Simulations studies of the di- 
amond melting curve have been reported for pressures up 
to 400 GPa [J, 1400 GPa and 2000 GPaQ. The last two 
studies were carried out by using "ab-initio" Molecular 
Dynamics. However, it would be prohibitively expensive 
to study nucleation using such an approach. We there- 
fore use a semi-empirical many-body potential that has 
been fit to experimentally measured and ab-initio calcu- 
lated properties of carbon solid phases and the liquid [f| . 
We use this model to study diamond nucleation in a sys- 
tem of 2744 particles in the "low-pressure" (PjlOOGPa) 
region of the phase diagram. In thispressure-range, the 
calculated melting line of Refs. [H, [ij are in reasonable 
agreement. In particular, all the three calculations pre- 
dict a melting temperature of about 7000 K at 100 GPa. 

In order to estimate the crystal-nucleation rate, we 
first determine the free-energy barrier AG cr it to form a 
critical nucleus. For state point A we use biased Monte 
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FIG. 1: The figure shows part of the carbon phase diagram 
from Ref. [I| and the iso-nucleation rate zones. The solid red 
lines represent the coexistence lines from. A is at Pa =85 GPa, 
T A = 5000 K, and B at P B = 30 GPa, T fl =3750 K. In the text, 
we give estimates for the nucleation rate at A and B. Along the 
green dashed curve, the ratio of 3-fold and 4-fold coordination 
in the liquid is 1:1. The color code used in the plot is: numbers 
on the right indicate the order of magnitude of the nucleation 
rate (in m _3 s _1 ). The continuous black curve indicates the 
boundary of the region where the nucleation rate is negligible 
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Carlo simulations ("umbrella sampling" [6J]) to estimate 
the excess Gibbs free energy of small crystallites, and 
find AG cr u a = 25 k B T for a critical nucleus size of 
Na = 110 pj. Knowing AG cr u and the kinetic pre- 
factor (see Ref. (H) we estimate the crystal nucleation 
rate to be Ra = 10 30 s _1 to -3 . When we tried to fol- 
low the same procedure to compute the diamond nucle- 
ation rate at state point 5(30 GPa, 3750 K), we failed 
to reach a critical nucleus that was small enough to fit 
in a system of 2744 particles. We therefore had to resort 
to an indirect way, based on CNT, to estimate AG cr i t . 
CNT assumes that AG (AT), the Gibbs free energy dif- 
ference between a metastable liquid containing an N- 
particle crystal nucleus and a pure liquid, is given by 
AG(N) = S(N)~f LS - N\A[i\, where S(N) is the area 
of the interface between an N-particle crystallite and the 
metastable liquid. In order to determine the number of 
particles in a crystallite, we use a spherical-harmonics 
based criterion (see e.g. Ref. Q) that allows us to distin- 
guish particles in a liquid-like environment from those in 
a crystal (diamond or graphite). The surface area S(N) 
is given by a(N/ ps) 2 ^ 3 , where the factor a depends on the 
geometry of the nucleus. From our simulations, we can 
only determine the product a^LS'- it is this quantity and 
the degree of supersaturation (Ap), that determine the 
nucleation rate. In order to calculate Ap, we compute 
the temperature dependence of the molar enthalpy differ- 
ence (Ah) between the supercooled fluid and the stable 



crystal at equal pressures. From Ah, Ap is evaluated by 
thermodynamic integration from the melting point [9(. 
We find: \Ap A /k B T\ = 0.60 and \Ap B /k B T\ = 0.77, re- 
spectively. 

^From the calculated AG cr it,A and the number density 
of the solid (pa= 0.191 A -3 ), we can estimate the surface 
free energy per unit area at state point A using Equa- 
tion([2|). Assuming that the critical nucleus is effectively 
spherical, we find 1ls ,a ~ 0.27 /c B T/A 2 =1.86 J/m 2 . We 
stress that, in what follows, we do not make use of this 
estimate: rather, we always employ the combination 07 
that follows directly from the simulations. 

In state point B we could not follow the same proce- 
dure, as a system of 2744 particles is too small to accom- 
modate a critical nucleus. In order to estimate Jls,b, we 
therefore prepared a rod-like crystal in a system with a 
slab geometry (a flattened box containing N ~ 4000 par- 
ticles, with lateral dimensions that are some four times 
larger than its height). The crystal rod is oriented per- 
pendicular to the plane of the slab. It spans the height 
of the simulation box and is continued periodically. The 
cross section of this crystal rod is lozenge shaped, such 
that its [lll]-faces are in contact with the liquid [lfj|. We 
used umbrella sampling to determine the Gibbs free en- 
ergy of such a crystallite as a function of its size, both 
at state points A and B. In this way we estimate the 
ratio of the surface free energies at A and B. We find 
that aj LS ,B/aiLS,A = Jls,b/jls,a ~ 2. 5 (the a factors 
are the same since the shape of the rod nucleus is ap- 
proximately the same at the two state points). Since we 
know 7ls,a from the height of the nucleation barrier in 
state point A for a spherical nucleus, we deduce the cor- 
responding 7ls,s f° r a spherical nucleus. Using our esti- 
mate, 7ls,a ~1.86 J/m 2 , we find ~/ ls ,b «0.68 k B T/k 2 = 
3.5 J/m 2 . As Afi and p B are known (p B =0.17 A~ 3 ), we 
can now use CNT to estimate AG cr ;t in state point B. It 
turns out that, mainly because Jls,b is 2.5 times larger 
than jls,a, the nucleation barrier in B is more than ten 
times higher than in point A, thereby hugely suppressing 
the nucleation rate (R B ~ 10~ 80 s _1 m -3 ). We can esti- 
mate the size of the critical (spherical) nucleus at B to 
be around N B =1700 particles. Thus, one would need a 
system of at least 17000 particles to contain a critical nu- 
cleus and avoid spurious interactions among its periodic 
images. Such a system size is beyond our present com- 
putational capacity. In contrast, in the slab geometry we 
find that the free energy of a lozenge-shaped crystal goes 
through a maximum at a size of ~ 340 particles, which 
is much less than the system size (4000 particles). 

To understand the microscopic origin for the large dif- 
ference in nucleation rates in state points A and B, it is 
useful to compare the local structure of the liquid phase 
in both state points. It turns out that the liquid struc- 
ture in state points A and B is markedly different (see 
also Refs. [ll|, 03) : liquid carbon is mainly four-fold 
coordinated at state point A (20% three-fold and 80% 
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FIG. 2: Typical snapshot of a small crystalline nucleus of 
~ 75 particles obtained at 30 GPa and 3750 K (B), sur- 
rounded by mainly three-fold coordinated liquid particles 
(gray lines). The nucleus contains both three- (left part) and 
four-fold (right part) coordinated particles. 



four- fold) , while at the lower temperatures and pressures 
of point B, the coordination in the liquid resembles that 
of the graphite and is mainly three-fold coordinated ( 5% 
two- fold, 85% three- fold and 10% four- fold ). Apparently, 
it is less favorable to create an interface between a dia- 
mond and a graphitic liquid than between a diamond and 
a four-fold coordinated liquid. The destabilizing effect of 
the graphitic liquid on the diamond nuclei is most pro- 
nounced for small nuclei (large surface-to- volume ratio). 
In fact, in state point B, nuclei containing less than 25 
particles tend to be graphitic in structure, with a small 
number of four-fold coordinated particles linking the dif- 
ferent graphite planes. Nuclei containing up to 60 parti- 
cles show a mixed graphite-diamond structure, whereas 
larger nuclei have a diamond bulk like structure, but the 
surface remains graphitic in nature (see Fig. 2). The 
unusual surface structure of the diamond nucleus is an 
indication of the poor match between a diamond lattice 
and a three-fold coordinated liquid. 

There are many network-forming liquids that, upon 
changing pressure and temperature, undergo profound 
structural changes or even liquid-liquid phase transi- 
tions |13l | - Our simulations on carbon indicate that such 
a change in the local coordination in the liquid has dra- 
matic consequences for the rate of crystal nucleation. Ex- 
periments on a completely different class of materials, viz. 
liquid metals [3], suggest that the local structure, in 
particular, local icosahedral packing, may interfere with 
direct nucleation of crystals. What is interesting about 
the present simulations is that we show that the ease of 
homogeneous crystal nucleation from one-and-the-same 
meta-stable liquid can be tuned by changing its pressure, 
and thereby its local structure. 

The thermodynamic conditions we discuss are relevant 
for experiments that study nucleation in compressed, 
laser-melted carbon. In addition, homogeneous nucle- 
ation of diamond may have taken place in carbon-rich 



white dwarf [15j. It has also been suggested that dia- 
monds could also have formed in the carbon-rich middle 
layer of Uranus and Neptune [l6| . The present work al- 
lows us to make a rough estimate of the conditions that 
are necessary to yield appreciable diamond nucleation on 
astronomical timescales. 

Neither white dwarfs nor planets consist of pure car- 
bon. Nevertheless, it is useful to estimate an upper 
bound to the diamond nucleation rate by considering 
the rate at which diamonds would form in a hypothet- 
ical environment of pure carbon. To this end we use 
our numerical data on the chemical potential of liquid 
carbon and diamond and our numerical estimate of the 
diamond-liquid surface free energy, to estimate the nu- 
cleation barrier of diamond as a function of temperature 
and pressure. We then use CNT to estimate the rate of 
diamond nucleation [l7| . The results are shown in Fig. 1. 
The figure shows that there is a region of some 1000 K 
below the freezing curve (continuous red line) where di- 
amond nucleation is less than 10 _40 m _3 s _1 . If the rate 
is lower than this number, not a single diamond could 
have nucleated in a Uranus-sized body during the life of 
the universe. As can be seen from the figure, our sim- 
ulations for state point B are outside the regime where 
observable nucleation would be expected. However, this 
figure provides just an upper bound to the rate of dia- 
mond nucleation. In practice, the carbon concentration 
is somewhat less in carbon-rich stars (~50%) [lH, and 
much less (1-2% \v$) in Uranus and Neptune. 

In Fig. 3, we show the effect of dilution on the regime 
where diamond nucleation is possible. To simplify this 
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FIG. 3: Diamond nucleation boundaries as a function of car- 
bon concentration: the rate is zero (no thermodynamic driv- 
ing force to nucleation) in the top region (liquid) , it is negligi- 
_1 ) in the middle region and non-negligible 



(> 10- 



in the bottom-right region. The nucleation 



rate is negligible when it corresponds to less than one nucleus 
per Uranus-sized planet over a period of 10 10 years. The 
left hand y-axis represents the temperature; the right-hand 
y-axis indicates the corresponding pressure for a Uranus-like 
isentrope (see Ref. [Tr|V 
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figure, we do not vary pressure and temperature inde- 
pendently but assume that they follow the adiabatic re- 
lation that is supposed to hold along the isentrope of 
Uranus [Hj]. We make the assumption that nucleation 
takes place from an ideal mixture [3] of C, N,0 and 
H [23| - Not surprisingly, Fig. 3 shows that dilution of 
the liquid decreases the driving force for crystallization 
to the extent that no diamond phase is expected for C 
concentrations of less than 8%. As before, there is a wide 
range of conditions where diamonds could form in princi- 
ple, but never will in practice. Assuming that, for a given 
pressure, the width of this region is the same as in the 
pure C case (almost certainly a serious underestimate), 
we arrive at the estimate in Fig. 3 of the region where 
nucleation is negligible (i.e. less than one diamond per 
planet per life-of-the- universe) . From this figure, we see 
that quite high Carbon concentrations (over 15%) are 
needed to get homogeneous diamond nucleation. Such 
conditions do exist in white dwarfs, but certainly not in 
Uranus or Neptune. Hence, our simulations indicate that 
it is extremely unlikely that diamonds could ever have nu- 
cleated from the carbon-rich middle layer of Uranus and 
Neptune. 
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